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We investigate the electronic density redistribution of rotated bilayer graphene under a perpendic- 
ular electric field, showing that the layers are actually coupled even for large angles. This layer-layer 
coupling is evidenced by the charge transfer on these structures as a function of the external voltage. 
We find an inhomogeneous excess charge distribution that is related to the moire patterns for small 
angles, but that persists for larger angles where the carriers' velocity is equal to that of single layer 
graphene. Our results show that rotated bilayer systems are coupled for all rotation angles. 



I. INTRODUCTION 

The electronic properties of graphene can be modi- 
fied by piling up a few layers, thus changing the behav- 
ior of the charge carriers due to the interlayer coupling. 
One of the main reasons behind the interest in bilayer 
graphene (BLG) is precisely that its low-energy prop- 
erties are different from those of monolayer graphene. 
Bilayer graphene with Bernal (or AB) stacking has mas- 
sive chiral fermions, and a gap can be opened by apply- 
ing an electric field. ^ In contradistinction, the dispersion 
relation of BLG with direct (or AA) stacking is linear, 
and remains gapless with an applied gate voltage. Even 
though Bernal stacking is the lowest energy arrangement 
for bulk graphite, in few-layer graphene other stack or- 
derings are possible: indeed, bilayer graphene with A A 
stacking has been observed,'^' as well as twisted bilayer 
graphene,^^''^ which consists of two adjacent graphene lay- 
ers where one of the layers is rotated with respect to the 
other, resulting in larger unit cells with intriguing elec- 
tronic properties.'^''^ 

As a matter of fact, there is an ongoing controversy on 
the electronic characteristics of rotated bilayer graphene. 
Several experiments on rotated few-layer graphene grown 
on SiC show an electronic behavior similar to that of 
single-layer graphene, with the same carriers' velocity 
as that of an isolated graphene monolayer; for this rea- 
son, these systems have b een c onsidered as composed of 
uncoupled graphene sheets .'^'^^Hllijjowever, experimental 
results in twisted graphene bilayers fabricated by chemi- 
cal vapor deposition indicate that twisted graphene bilay- 
ers may be coupled, especially for small rotation angles, 
as evidenced by the appearance of low-energy van Hove 
singularitie^i^ and the measured renormalization of the 
carrier velocity.li^ 

From the theoretical viewpoint, it is established that 
twisted BLG with a relative rotation angle (RRA) greater 
than 10° presents a linear dispersio n rela tion with the 
same velocity as monolayer graphene.'^'^'^ For RRA be- 
tween 1° and 10°, the carrier velocity diminishes, as evi- 
denced by continuunP and combined tight-binding/ first- 



principles calculations.'^' For smaller angles, around 1°, 
flat bands appear close to the Fermi energy.^ ^ However, 
for certain rotation angles where theoretical calculations 
predict a renormalization of the carriers' velocity, some 
experimental results hav e shown a behavior like that of 
monolayer graphend^iEl while others do obtain the pre- 
dicted renormalized velocity for small RRA.ES In princi- 
ple, the key to understand this disagreement is to eluci- 
date the coupling strength in adjacent graphene layers.l^ 

In this paper we propose a way to assess the interlayer 
interaction in twisted bilayer graphene by exploring the 
spatial distribution of the electronic density. We show 
that, in spite of the linear dispersion relation and the 
Fermi velocity of these systems, the layers can be actually 
coupled even for large rotation angles. The layer-layer 
coupling is evidenced by the charge transfer between lay- 
ers under an applied electric field. In addition, we find an 
inhomogeneous excess charge for small angles, persisting 
for larger angles where the carriers' velocity is equal to 
that of single layer graphene. We have found that the av- 
erage excess charge density on a layer of a twisted bilayer 
graphene, due to an applied electric field, depends on the 
RRA. The charge transfer is larger the smaller the angle 
but even for RRA above 20° there is a significant charge 
transfer. The charge has an inhomogeneous distribution 
over the twisted unit cell, with a relative weight on sites 
A and B also depending on the angle. 

This article is organized as follows. In Section |TT] we 
describe the geometry of the commensurate graphene bi- 
layers and the model employed to calculate the electronic 
properties. Section [lll| is devoted to explain and discuss 
our results. We conclude in Sec. |IV[ where we summarize 
our main findings. 



II. GEOMETRY AND MODEL 

A. Geometry 

The system studied consists on two graphene layers 
which are rotated from a Bernal stacking configuration. 



2 



subject to a potential difference. We consider the top 
layer to be at a positive potential +V and the bottom 
layer at —V. This potential difference produces an elec- 
tronic charge transfer from the bottom to the top layer. 
Consequently, an electronic excess density n builds up 
on the top layer. We will show that the excess charge 
depends on the interlayer coupling. 

In order to build a commensurate unit cell we follow 
a procedure similar to that of Campanera et al, find- 
ing coincidence lattice points in the BLG crystal.^- We 
start with two graphene layers with Bernal (AB) stack- 
ing. For this arrangement there are two non-equivalent 
sites within a layer, namely, the A site, where an atom 
of the upper layer lies directly on top of another atom of 
the lower layer, and the B site, where the atom is exactly 
at the center of the hexagon of the lower layer. We have 
chosen a B site as our rotation center. We do a clockwise 
commensurable rotation from a vector r — mdi + na2 
to ti — nai + 1710,2, where ai — (— 1/2, •\/3/2)ao and 
0.2 = (1/2, a/3/2) flo are the graphene bilayer lattice vec- 
tors; n,m are integers, and oq = 2.46 A is the lat- 
tice constant. The unit cell vectors can be chosen as 
ti — ndi + 7X10,2 E^nd ^2 = — mSi + in + m)a2. All 
magnitudes needed, like the relative rotational angle 
(RRA), the number of atoms in the unit cell and the 
reciprocal lattice vectors, can be expressed as functions 
of 71, m and oq. For instance, the RRA is given by 
cosO = (n^ + Amn + m^)/2(7i^ -|- 7nn + m^). Hence- 
forth, we label a commensurate unit cell by the indices 
(n, m), which determine uniquely the atomic structure of 
the corresponding rotated bilayer. 

Fig. [l] shows the unit cell of the (4, 3) bilayer, along 
with the unit cell vectors ti and 12- In a rotated bilayer, 
we distinguish four distinct regions, centered in specific 
sites given by the different possible atomic stackings. We 
label as AB a site with a top atom at the center of an 
hexagon. An A A site has an atom exactly on top of 
another, with a similar stacking for the nearest neighbors, 
and a BA site has an atom in the bottom layer exactly 
at the center of an hexagon of the top layer. These sites 
are indicated in Fig. [T] 

The unit cells presenting distinct moire patterns cor- 
respond to small angles (below 10°); therefore, in order 
to obtain large non-trivial unit cells, we choose n and 
m as two integers without any common divisors. By 
varying these numbers, any RRA from 0° to 60° can 
be attained.Ii^l It should be noted that structures with 
rather close RRA can be obtained from very different n 
and TJi. Such structures have similar moire patterns, but 
they arise from unit cells with quite different number of 
atoms. We did most of our calculations for m — n — 1. 
This choice gives the smallest unit cell for a set of RRA 
and covers fairly well the range of angles of interest . How- 
ever, other structures with dissimilar n, m indices were 
generated in order to check that our results do not de- 
pend on the shape of the unit cell. 

For a bilayer graphene or graphite, moire patterns are 
often observed for angles below 10°. For instance, a 




FIG. 1: A twisted unit cell (4,3) for a RRA of 9.43°; it con- 
tains 148 atoms. The AA, AB and BA sites are shown. 

(17,16), with a 2° RRA, results in a unit cell with 3268 
atoms. Setting the rotation axis on a B site and a unit 
cell with TO = n — 1 gives a unique AA site and an AB 
site along the diagonal of the unit cell. The AA site lies 
at 1/3 of the diagonal and the AB site lies at 2/3. From 
our construction we get at every corner of our unit cell 
a BA site. There are also some other interesting points 
inside the unit cell; we will call a slip site (SL)^^ a point 
like that in the middle of the line joining an AB site and 
the end of the main diagonal line. A slip site apparently 
looks as derived by a relative translation of the two lay- 
ers, hence its name. 



B. Tight Binding Model 

We model the bilayer graphene band structure within 
the tight binding approximation. Within each layer, we 
consider a fixed nearest-neighbor intralayer interaction 
7o — 3.16 eV. For the layer-layer interaction we consider 
a distance-dependent hopping between more than near- 
est layer-to-layer neighbors. Thus, the hamiltonian is 
given by H = Hi + H2 + Hint , where Hi and H2 are the 
hamiltonians for the top and the bottom layer and Hint 
describes the interlayer coupling, 

H.nt = -J2 lie-^^'"''-''^clc, + h.c, (1) 

where 71 = 0.39 eV is the nearest-neighbor interlayer 
hopping parameter, d is the interlayer distance, rjj is the 
distance between atom i on the top layer and atom j on 
the bottom layer, and /3 = 3. This value of /? reproduces 
accurately the dispersion bands calculat ed with a Den- 
sity Functional Theory (DFT) approach! ^ ' ^^ ' ^^ ' Besides, 
it gives a value of 73, the second nearest neighbor in- 
terlayer hopp ing, in agreement with values employed by 
other authors .E^Efil We let every atom in the top layer to 
interact with the atoms in the bottom layer located inside 
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a circle of radius Sao, taking into account the complex- 
ity of the unit cell. This breaks the electron-hole (e-h) 
symmetry due to the fact that we are mixing the two sub- 
lattices. The asymmetry depends on the value of /3 and 
on the RRA, being larger for smaller angles; it manifests 
in the difference between the electron and hole velocities. 
We find a larger velocity for electrons, as in the experi- 
mental results obtained by Luican et al.^^with a greater 
e-h asymmetry in the experimental values. 

In order to compute the charge density and total en- 
ergy, we perform a sum over the Brillouin Zone (BZ); the 
k points were selected with a Monkhorst-Pack scheme,!^ 
taking into account the symmetries of the BZ. Conver- 
gency tests were made to select the appropriate number 
of k points for the mesh. 

Under an external bias, the LDOSJ^ and the charge 
density in a twisted unit cell are inhomogeneous. This, 
along with the discreteness of the crystal lattice points, 
lead us to compute the Fourier transform (FT) of the 
charge excess density. As there are large oscillations of 
the charge density between the two sublattices, we com- 
puted separately the FT contributions of the A and B 
atoms. Furthermore, the FT allows us to evaluate in a 
continuous way the charge density on AB, BA, and A A 
regions on the top layer. Only the zero order and the 
three lower components were needed to fully reproduce 
the charge distribution. 



III. RESULTS AND DISCUSSION 

As discussed in the preceding section, when a poten- 
tial difference is applied to a coupled bilayer, charge is 
transferred from one layer to the other. In the geometry 
chosen, there is an excess electronic charge n that builds 
up on the top layer. Note that if the layers are uncou- 
pled, there is also a charge transfer due to the relative 
shift of the Dirac cones in the two layers, but the amount 
will be smaller than in a bilayer with nonzero coupling. 

In Fig. [2] we plot the excess electronic density on the 
top layer as a function of V for several rotated bilay- 
ers, including the AB case for comparison. We see how 
the excess electronic density diminishes with increasing 
rotation angle, indicating that the coupling between the 
rotated graphene layers is reduced for higher RRA, albeit 
non-negligible in any case. The values for large angles are 
rather similar; this saturation is correlated to the behav- 
ior of the carriers' velocity for increasing angle, which is 
equal to that of monolayer graphene for large RRA, see 
inset of Fig. |2] 

In Fig. [3] we show the angle dependence of the charge 
transfer for a fixed, nonzero electric field {V = 0.08 V). 
The excess electronic charge has a strong variation for 
small angles, and then tends to a constant value, al- 
though there is a smooth variation even in the range from 
20° to 30°. However, this limiting value for the excess 
electronic charge is much higher than the one obtained 
for an uncoupled bilayer: the inset of Figure |3]depicts the 
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FIG. 2: (Color online) Excess electronic density on the top 
layer as a function of V for several rotated bilayer graphene 
structures. The rotation angles are indicated in the Figure. 
Inset: Fermi velocity in twisted bilayer graphene, vt, nor- 
malized with respect to that of monolayer graphene, vq, as a 
function of the rotation angle RRA. 

excess charge on a bilayer AB subject to the same elec- 
tric field as a function of the interlayer coupling constant. 
We label this fictitious varying interlayer hopping 71 to 
distinguish it from the parameter 71 defined in Section 
|IIB| We see that the charge transferred between layers 
in an uncoupled system (71 — 0) is n ^ 0.21 x 10^^ cm^^, 
roughly 30% smaller than the excess electronic charge for 
the largest RRA value in Fig. [3] almost n = 0.3 x 10^^ 
cm^^. Therefore, reaching a saturation in the charge 
transfer does not necessarily mean that the twisted bi- 
layers are uncoupled. Notice as well that the charge with 
diminishing angle tends to the value of the perfect stacked 
AB bilayer. 

The excess charge presents an inhomogeneous spatial 
distribution that depends on the rotation angle. As there 
is a strong oscillation in the density between neighbor- 
ing A and B sites, we choose to plot it for these points 
separately. In Fig. [4] we show the excess density on the 
top layer for a (15,14) bilayer. Different regions with re- 
spect to the electronic density distribution can be distin- 
guished, corresponding to the stacking regions previously 
identified (see Fig. [l]). For example, the BA zones, lo- 
cated at the corners of the unit cell, present an increase 
in the electronic density in the A atoms, whereas it de- 
creases on the B atoms. On the other hand, AB zones 
display the opposite behavior: the A atoms show a re- 
duction of the electronic charge, which in turn increases 
on the B atoms. The asymmetry is quite large, as can 
be seen in Fig. l4] The excess charge density in B (A) 
atoms in AB (BA) sites is 300 % larger than the av- 
eraged value (n — 0.7 x 10^^ cm^^, see Fig. [3|, varying 
from -1.4 to 2.8 x 10^^ cm"^. On the other hand, the AA 
regions present smaller charge fiuctuations, with similar 
electronic density values for A and B atoms. This can be 
understood by considering the local surroundings of the 
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FIG. 3: (Color online) Excess electronic density on the top 
layer as a function of the rotation angle for a fixed electric 
field {V = 0.08 V). Inset: Excess charge density as a function 
of the interlayer coupling for AB bilayer graphene. 



atoms: the AB or BA regions resemble the stacked AB 
bilayer, where there is a difference between the LDOS 
and the charge density in both sites, whereas the A A 
region is more like the stacked AA bilayer, with similar 
density values in A and B sites. 
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FIG. 4: (Color online) Contour plot of the excess electronic 
density, in units of 10^^ cm~^, on the top layer in the A (left) 
and B (right) atoms for the (15, 14) rotated bilayer, with 2524 
atoms in the unit cell. The position coordinates {x and y axes) 



Given such differences, we have studied the charge dis- 
tribution in different stacking regions as a function of the 
rotation angle in the top layer. We average the FT of the 
excess charge density in atoms A and B, evaluated at the 
AA, AB, BA, and shp points of the unit cell. Fig. [5] 
shows this averaged FT for a fixed V — 0.08 V at the 
AA, AB and slip regions for angles below 10°. As there 
is a symmetry under the interchange of A and B sites 
on the AB and BA regions, only the values of the AB 
site are shown. Several features should be noted: (i) the 
AA site takes less charge than the AB site; (ii) the slip 



site density behavior is very similar to that of the total 
density; and (iii) for very low angles, the density is more 
homogeneous, as opposed to the LDOS, which peaks at 
AA sites near the Fermi energy for small RRA.^^ 
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FIG. 5: (Color online) Excess electronic density on the top 
layer as a function of the rotation angle for a fixed potential 
V — 0.08 V for different stacking regions on the unit cell, 
namely, AA, AB and slip. 

Besides these magnitudes related to the charge trans- 
fer under an applied perpendicular electric field, one way 
to theoretically assess the electronic coupling is to cal- 
culate the interlayer coupling kinetic energy in the un- 
biased double layer. We define Ecoupi as the difference 
between the total kinetic energy of the coupled bilayer 
and the kinetic energy of the uncoupled double layer, 
Ecoupi = {'^\Hint\'4^) , where j-f/)) is the ground state wave- 
function of the total coupled bilayer hamiltonian. This 
quantity is directly related to the interlayer interaction, 
which should vary with the RRA. Fig. [6] presents the 
interlayer coupling energy Ecoupi as a function of the ro- 
tation angle. We see that although Ecoupi increases with 
increasing rotation angle, is by no means negligible for 
RRA over 20°. The energy for low angles decreases, tend- 
ing to the value of the stacked AB bilayer, —5.1 meV. 

Finally, we note that in the analysis of the bias-induced 
charge transfer we have neglected the Co ulomb energy 
Ec corresponding to the two charged layers .1221231 We can 
estimate the charging energy Ec by considering that the 
bilayer is like a capacitor with oppositely charged plates, 
and assuming that the electronic density is uniformly 
distributed in the two-dimensional layers. With these 



assumptions, Ec — ^-^^do, where e is the elementary 
charge, dg is the interlayer separation, S is the sample 
area, and e is the dielectric constant of the medium. This 
Coulomb energy induces an extra potential between the 
layers given by A]/(n) = . Because of the small sep- 
aration between the two graphene layers, the Coulomb 
potential induced by the charge transfer is much smaller 
than the externally applied bias, so therefore it is justi- 
fied to neglect it in the calculation of the bias-induced 
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charge imbalance. 
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means. Firstly, the application of an external electric 
field induces a charge transfer between layers that can 
be explored as a way to indicate the electronic coupling, 
this charge transfer depends on the RRA. Secondly, the 
spatial charge inhomogeneities are also a signal of the 
electronic coupling between layers, analogously to the dif- 
ferences in the LDOS reported by other authors.'^^ The 
coupling was theoretically estimated by computing the 
coupling energy, i.e., the difference in the energy of a cou- 
pled bilayer and the energy of two uncoupled graphene 
layers. This magnitude is nonzero for all angles, thus ev- 
idencing that rotated graphene bilayers are coupled for 
all rotation angles. The values of the excess charge den- 
sity and the coupling energy tends toward the stacked 
AB bilayer for low twist angles. 



FIG. 6: (Color online) Total coupling energy per atom of the 
twisted bilayers as a function of the RRA. 



Acknowledgments 



IV. CONCLUSIONS 

In summary, we have shown that the coupling between 
rotated graphene bilayers can be evidenced by several 



This work has been paxtially supported by MEC-Spain 
under grant FIS2009-08744 and by the CSIC/CONICYT 
program, grant 2009CL0054. E.S.M. acknowledges CON- 
ICYT (Chile) and UTFSM for the internal grant PIC- 
DGIR P.V. acknowledges FONDECYT grant 1100508. 



^ E. McCann and V. I. Fal'ko, Phys. Rev. Lett. 96, 086805 
(2006). 

^ Z. Liu, K. Suenaga, P. J. F. Harris, and S. lijima, Phys. 
Rev. Lett. 102, 015501 (2009). 

^ J. Hass, F. Varchon, J. E. Millan-Otoya, M. Sprinkle, 
N. Sharma, W. A. de Heer, C. Berger, P. N. First, L. Ma- 
gaud, and E. H. Conrad, Phys. Rev. Lett. 100, 125504 
(2008). 

* T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Roten- 
berg. Science 313, 951 (2006). 

^ A. Reina, X. Jia, J. Ho, D. Nezich, H. Son, V. Bulovic, 
M. S. Dresselhaus, and J. Kong, Nano Letters 9, 30 (2009). 

® J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. 
Castro Neto, Phys. Rev. Lett. 99, 256802 (2007). 
E. Suarez Morell, J. D. Correa, P. Vargas, M. Pacheco, and 
Z. Barticevic, Phys. Rev. B 82, 121407 (2010). 

* G. Trambly de Laissardiere, D. Mayou, and L. Magaud, 
Nano Letters 10, 804 (2010). 

^ R. Bistritzer and A. H. MacDonald, PNAS (2011). 

° F. Varchon, P. Mallet, L. Magaud, and J.-Y. Veuillen, 

Phys. Rev. B 77, 165415 (2008). 
^ Y. J. Song, A. F. Otte, Y. Kuk, Y. Hu, D. B. Torrance, 

P. N. First, W. A. de Heer, H. Min, S. Adam, M. D. Stiles, 

et al., Nature 467, 185 (2010). 
^ J. Hicks, M. Sprinkle, K. Shepperd, F. Wang, A. Tejeda, 



A. Taleb-Ibrahimi, F. Bertran, P. Le Fevre, W. A. de Heer, 
C. Berger, et al., Phys. Rev. B 83, 205403 (2011). 
L. Guohong, A. Luican, J. M. B. Lopes dos Santos, A. H. 
Castro Neto, A. Reina, J. Kong, and E. Y. Andrei, Nature 
Physics 6, 109 (2010). 

A. Luican, G. Li, A. Reina, J. Kong, R. R. Nair, K. S. 
Novoselov, A. K. Geim, and E. Y. Andrei, Phys. Rev. Lett. 
106, 126802 (2011). 

S. Shallcross, S. Sharma, and O. A. Pankratov, Phys. Rev. 
Lett. 101, 056803 (2008). 

A. H. MacDonald and R. Bistritzer, Nature 474, 453 
(2011). 

J. M. Campanera, G. Savini, L Suarez-Martinez, and M. I. 
Heggie, Phys. Rev. B 75, 235449 (2007). 
S. Latil, V. Meunier, and L. Henrard, Phys. Rev. B 76, 
201402(R) (2007). 

B. Partoens and F. M. Peeters, Phys. Rev. B 74, 075404 
(2006). 

^° D. Chung, Journal of Materials Science 37, 1475 (2002). 
H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 
(1976). 

^2 E. McCann, Phys. Rev. B 74, 161403 (2006). 

J. Fernandez-Rossier, J. J. Palacios, and L. Brey, Phys. 
Rev. B 75, 205441 (2007). 



